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Abstract 

We report the development and implementation 
of hybrid methods that combine Quantum Me- 
chanics (QM) with Molecular Mechanics (MM) to 
theoretically characterize thiolated gold clusters. 
We use, as training systems, structures such as 
Au 25 (SCH 2 -R)i8 and Au3 8 (SCH 2 -R)24, which can 
be readily compared with recent crystallographic 
data. We envision that such an approach will lead 
to an accurate description of key structural and 
electronic signatures at a fraction of the cost of a 
full quantum chemical treatment. As an example, 
we demonstrate that calculations of the 1 H and 
13 C NMR shielding constants with our proposed 
QM/MM model maintain the qualitative features 
of a full DFT calculation, with an order-of-magnitudc 
increase in computational efficiency. 

1 Introduction 

It is well known that gold nanoparticles are capa- 
ble of catalyzing a number of fundamental chemi- 
cal reactions [TJ HJ |3] . Tunable reactions in metallic 
nanoparticles can be obtained by encapsulation of 
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nanocrystals in thiolated monolayers [H [5] known 
as monolayer protected clusters (MPCs). Exploit- 
ing MPCs for applications in the area of catalysis 
requires a proper structural and dynamical char- 
acterization of the physical properties of the pro- 
tecting layer. It was not until recently, that such 
requisite has been met by the total structure deter- 
mination of a number of thiolated gold clusters [51 
[7J IB] • Such structures are characterized by surface 
gold atoms in the so called "staple" (-S-Au-S-) or 
"V-shape" motifs (-S-Au-S-Au-S-) , as opposed 
to commonly assumed structures in which thio- 
lates only passivate a high symmetry gold cluster. 
Along these experimental studies, Density Func- 
tional Theory has been crucial to interpret struc- 
tural data under the light of NMR spectroscopy 
[HI HO] ; optical absorption [TT] , and electrochemical 
experiments [12]. Due to the size of these clusters, 
augmented by thiolated ligands, theoretical char- 
acterization (e.g calculation of spectra) requires 
an enormous amount of computational resources 
(both in time and memory). Intrinsic limitations of 
current Quantum Chemistry software places these 
types of calculations at the limits of computational 
tractability. This is clearly the case if, in addition 
to quantum detail, finite temperature simulations 
are required. In fact, a number of studies have 
shown a marked dependence of optical and elec- 
trochemical properties on temperature [T31 [2] [T5] . 
Thus, to retain quantum detail and incorporate fi- 
nite temperature motions it would certainly be ad- 
vantageous to develop approaches that can com- 
bine these features with a reasonable computa- 
tional cost. 



Empirical potentials have been used in the past 
for gold clusters to search structural minima fol- 
lowed by quantum mechanical refinement |16) , but 
little has been done in the context of MPC's, es- 
pecially considering the recent re-evaluation of the 
sulfur-gold binding motifs. Another approach would 
be to develop hybrid methods such as QM/MM, 
which were originally developed in the context of 
enzyme reactions [T7] . In this work we develop MM 
and QM/MM models, explore different implemen- 
tations, and evaluate their ability to reproduce ex- 
perimental data. We envision that such approaches 
will be of value to enable, in a practical manner, 
sampling of a large number of monolayer confor- 
mations. In particular, such procedures will be rel- 
evant to describe highly interacting oligopeptide 
ligands [TSUH]- In addition, partitioning the system 
into a QM and MM region will lead to several ad- 
vantages. Apart from the obvious computational 
benefits, such method opens up a wide range of 
possibilities for certain studies of similar systems 
for which a full QM description is impractical and 
often unnecessary This new approach also opens 
the possibility of focusing strictly on the electronic 
structure of the monolayer alone (or of the gold 
cluster and a fraction of the monolayer) which, 
for instance, can be advantageous to obtain insight 
into catalytic properties of these MPC's[TH|. Thus, 
the goal of this study is to develop a first gener- 
ation of hybrid methods that can both decrease 
the computational cost and, at the same time, re- 
produce the intricacies of the various gold-sulfur 
surface motifs. 

We use [Au25(SCH3)i 8 ]~ as the central proto- 
type model and extend our theory to explain struc- 
tural and chemical properties of larger MPC's. We 
begin by determining a minimum set of parameters 
and functional forms that would be sufficient to de- 
scribe a thiol-protected gold cluster via a Molecu- 
lar Mechanics (MM) force field. We then use these 
parameters to develop a hybrid QM/MM model 
which treats regions of chemical interest at a DFT 
level of theory, keeping all other atoms in the clus- 
ter at the MM level. 



2 Computational Methods 

We propose the use of hybrid QM/MM models to 
study physical and chemical properties of mono- 
layer protected gold clusters. Au-SCH 3 , Au-(SCH 3 )- 
Au, and [Au25(SCH 3 )i 8 ]~ are used as the training 



systems. In the following sections, we present the 
parameters and functions necessary to describe the 
structure in a molecular mechanics (MM) frame- 
work. We continue our discussion by introducing 
QM/MM hybrid models to predict or refine struc- 
tures and NMR properties of similar gold clusters. 
All QM/MM calculations were performed with the 
Gaussian 09 suite [5^ using the two-layer ONIOM 
scheme [3]] [55]. -"- n this scheme, the entire system is 
divided into two regions (X and Y). The QM/MM 
energy is obtained via an extrapolation of three in- 
dependent calculations: 

E = E(QM) X + E(MM) X +Y - E{MM) X (1) 

where E(QM)x is the energy of region X at the 
QM level, E(MM)x+y is the energy of the entire 
system (X and Y) at the MM level, and E(MM) X 
is the energy of X at the MM level. In the so called 
"electronic embedding" (EE) approach, electrostatic 
interactions between X and Y are included in each 
of the terms of the right hand side of Eq. [l] so 
that electrostatic interactions are canceled out at 
the MM level, but remain at the QM level. In the 
"mechanical embedding" (ME) approach, electro- 
static interactions between region X and Y are 
only included at the MM level (the last two terms 
of Eq. [l). In both EE and ME, Van dcr Waals 
(VDW) interactions between X and Y are consid- 
ered at the MM level (in the term E{MM) x +y)- 
All calculations presented here make use of the EE 
approximation. 



2.1 Molecular Mechanics Model 

Construction of QM/MM models was made un- 
der the assumption that even a pure MM force 
field description of the thiolated clusters should re- 
produce the experimental structural data, at least 
in a qualitative and semi-quantitative manner. To 
this end, we sought to develop a force field with 
a minimal set of parameters and functional forms. 
Derivation of the force field was guided by the no- 
tion that these clusters present two well defined 
domains. Domain 1 contains the inner core (i.e. 13 
Au atoms for Au25(SCH2-R)is and 23 Au atoms 
for Au38(SCH 2 -R)24) forming high symmetry packed 
structures. For this domain we defined a 6-12 Lcnnard- 
Jones potential as implemented in the Amber force 
field [23 . Domain 2 contains the thiolated ligands 
and all Au atoms involved in the V-shape or sta- 
ple motifs (Fig. [IJ. All gold atoms in this domain 
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Fig. 1 [COLOR ONLINE] Possible binding patterns 
between gold and thiolated ligands 



Table 1 6-12 Lennard Jones parameters for a gen- 
eral MPC. Partial charges correspond to thiol in 
[AuastSCHajis]- 1 



Atom 


Atom-type 


cr 


e 


charge 


Au 


Au 


1.5000 


1.0000 


0.00000 


Au 


AuC 


1.5000 


1.0000 


0.00000 


Au 


AuS 


1.9000 


1.0000 


0.00000 


S 


SS 


2.0000 


0.2500 


-0.22417 


S 


SC 


2.0000 


0.2500 


-0.22417 


C 


CT 


1.9080 


0.1094 


-0.15441 


H 


HI 


1.3870 


0.0157 


0.10767 



Table 2 Stretching parameters corresponding to the 
force field formula (l/2)fc e (r — r eq ) . 

At-typel At-type2 k r r eq 



AuC 


SC 


150.00 2.48 


AuS 


SC 


150.00 2.38 


AuS 


SS 


150.00 2.38 


SC 


CT 


237.00 1.84 


SS 


CT 


237.00 1.84 


CT 


HI 


340.00 1.09 


CA 


HI 


340.00 1.09 



Table 3 Bending parameters corresponding to the 



force field formula (l/2)k$ (0 - 



i ) 



At-typel 


At-type2 At-type3 


kg 




AuC 


SC CT 


35.00 


105.00 


AuS 


SC CT 


35.00 


105.00 


SC 


CT HI 


50.00 


109.50 


HI 


CT HI 


35.00 


109.50 


SC 


AuS SS 


20.00 


180.00 


AuS 


SS CT 


35.00 


105.00 


AuS 


SS AuS 


20.00 


100.00 


AuC 


SC AuS 


20.00 


91.00 


SS 


CT HI 


50.00 


109.50 


SC 


CT CT 


50.00 


111.50 


SS 


CT CT 


50.00 


111.50 


HI 


CT CA 


35.00 


110.50 


CA 


CA HI 


35.00 


120.00 



are assigned the same force field atom type (AuS). 
As starting values we took those reported in the 
UFF force field of Rappe et al [24] . These param- 
eters were modified systematically and stepwise 
until the geometry of [Au25(SCH 3 )i 8 ] _1 was re- 
produced within a tolerance of 0.2 A in the root 
mean square deviation. Gold atoms in 1 are sepa- 
rated from those that are only bonded to other core 
atoms (type Au) and those that are connected to 
core atoms and sulfur atoms (type AuC) . Although 
the VDW parameters are identical for these two 
types, this distinction is required for the definition 
of bonded parameters (vide infra). VDW parame- 
ters are reported in Table [l] 

Electrostatic interactions between the ligands 
and gold atoms in domain 1 are assumed to be neg- 
ligible. Support for this assumption is given by a 
Bader charge analysis [23 H5] of [Au 25 (SCH 3 ) 18 ]- 1 
which shows that charges in domain 1 are very 
small (« — 0.04e for the central core atom and 
rs 0.02e for the other core atoms), adding up to 
a total charge of 0.15e. Bader charges for AuS 
atoms are also small (sa 0.04e). Thus, to facili- 



tate transferability among different core sizes, the 
charges on all gold atoms are assumed to be zero. 
Charges on the thiol ligands are taken from an 
Electrostatic Potential (ESP) analysis (as imple- 
mented in Gaussian 09) for the neutral molecule 
Au2-[SCH3]2- This analysis provides the follow- 
ing charges: qs — —0.23605, qc = —0.16628, and 
<7h = 0.09580. To ensure that the total MM charge 
adds up to the total charge of the cluster (Q), a 
simple correction Sq is applied to each charge via 
the formula Sq = (Q/M — (qs + qc + 3<lH))/5, where 
M is the number of thiolated ligands. Table [T] re- 
ports these charges for Q = — 1 and M = 18, which 
uses the correction Sq = 0.01188. 

Bonded parameters (stretching, bending, tor- 
sion) involving X-Au-S-X were derived by fitting 
the force constants in all functional forms to match 
vibrational frequencies of selected modes of vibra- 
tion for the molecules Au-SCH 3 and (Au-SCH 3 - 
Au). All other parameters intrinsic to SCH 3 were 
taken from the Amber force field without modifi- 
cation. Tables[2][3j and[4]report all bonded param- 
eters involving the newly defined atom types and 
the existing types in the Amber force field. 
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Table 4 Amber torsional parameters corresponding to the force field formula 2Zi=i Mi [1 + cos(id — Oi(i + 4))]/N p , 



A-i- funo 1 






/vc-ty pe^t 






Ci 
<-J3 


U4 


A f 




A if 


A 1 
ill 4 


iVp 


AuC 


SC 


CT 


HI 














0.00 


0.00 


0.16 


0.00 


1.00 


AuC 


SC 


CA 


HI 





{) 








0.00 


0.00 


0.16 


0.00 


1.00 


AuS 


SC 


CT 


HI 














0.00 


0.00 


0.16 


0.00 


1.00 


AuS 


SC 


CA 


HI 














0.00 


0.00 


0.16 


0.00 


1.00 


AuS 


SS 


CT 


HI 














0.00 


0.00 


0.16 


0.00 


1.00 


AuS 


ss 


CA 


HI 














0.00 


0.00 


0.16 


0.00 


1.00 


SC 


AuS 


SC 


CT 














0.00 


0.00 


0.00 


0.00 


1.00 


SC 


AuS 


SC 


CA 














0.00 


0.00 


0.00 


0.00 


1.00 


SC 


AuS 


SS 


CT 
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0.00 


0.00 


0.00 


1.00 


SC 


AuS 
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0.00 
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0.00 


0.00 
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0.00 


0.00 


0.00 


0.00 


1.00 
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0.00 


0.00 


0.00 


1.00 




Fig. 2 [COLOR ONLINE] Structure of 
[Au25(SCH3)i8] - , with Au atoms shown in or- 
ange, S in yellow, C in gray. H atoms are not shown 
for clarity. The AU13 core (MM region) is shown with 
spheres while all other atoms (QM region) are shown 
with sticks. 



2.2 QM/MM Models 
2.2.1 Au2 5 

One of the major objectives in developing a full 
MM force field is to use it as the underlying po- 
tential in hybrid QM/MM models. We benchmark 
our model by comparing structural descriptors in 
[Au 2 5(SCH 3 )i 8 ]-. The crystal structure of [N(C 8 H 17 ) 
] + [Au25(SCH 2 CH 2 Ph) 18 ]~ 1 was determined recently 
by Heaven et al. [7] , which was followed by several 
other experimental investigations |27j and theoret- 
ical [5S] calculations using density functional the- 
ory (DFT) methods. The initial geometry for our 
QM/MM calculations is derived from that crystal 



structure by replacing the ethylphenyl groups with 
methyl groups. We adjust the resulting [Au25(SCH 3 )i 8 ]~ 
structure to obtain a symmetric structure belong- 
ing to the point group C,. 

The [Au25(SCH 3 )i 8 ]~ structure consists of an 
icosahedral Aui 3 core, protected by six V-shape 
motifs in an approximate octahedral arrangement. 
Each of the S atoms are connected to the organic 
ligand (CH 3 ), forming a monolayer protected clus- 
ter. We propose a QM/MM model in which the 
icosahedral Aui 3 core is treated via MM whereas 
the V-shape motifs with the ligands are treated via 
QM methods (see figure[2|. These two regions inter- 
act with each other via the ONIOM extrapolation 
scheme with electronic embedding as described in 

Eq.rjrj 

Our choice of the partition for [Au25(SCH 3 )i 8 ]~ 
cluster was fairly straight-forward, making use of 
the intuitive definition of domains presented in 
the previous section. The icosahedral Aui 3 core is 
very stable and can be accurately described by a 
Lennard- Jones potential. Each of the gold atoms 
in the core, except the central one, is bonded to a 
sulphur atom forming a V-shape motif. A "cut" is 
defined between atoms types AuC and SC. These 
bond cuts are implemented in ONIOM by defining 
a "link-atom" to satisfy the valence and avoid bro- 
ken bonds. The default choice of such a link-atom 
is typically a hydrogen atom. However for our pur- 
poses, the best choice of link-atom turned out to be 
another Au atom, thus forming Au(link)-S bonds 
in place of Au(core)-S bonds. Since the core has a 
small charge (ss 0.15e), and since the QM region 
needs to be assigned an integer charge, we assume 
a -1 charge solely distributed among the atoms 
of domain 2 (QM region), thus making the entire 
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Fig. 3 [COLOR ONLINE] Structure of Au 3 s(SCH 3 )24. 
The AU23 core (MM region) is shown with spheres while 
all other atoms (QM region) are shown with sticks. 



Table 6 13 C NMR chemical shifts S (in ppm) for 
[Au 2 5(SCH3)i2(SCH2CH 2 Ph) 6 ]- cluster. Reported are 
the calculated values within the QM/MM approxima- 
tion, full DFT calculation, and experimental values. 



Atom 


QM/MM 


DFT 


Experiment |10| 


<^ C inner 


50.75 


55.80 




Pouter 


A Q 71 


A(\ *7 A 


(17 

oo.y / 


ft C in ner 


45.11 


44.26 




ft C outer 


44.07 


43.47 


42.18 


i C inner 


137.35 


141.05 


141.88 


1 C outer 


137.43 


140.42 


140.88 


1X1— C in ner 


125.48 


125.00 


128.52 


FT1 Cottier 


125.61 


124.94 


128.63 


O Conner 


125.29 


126.49 


129.85 


O Co liter 


125.37 


125.42 


129.30 


P Ci 7i ner 


123.99 


122.45 


126.27 


P C outer 


123.96 


122.35 


126.49 



MUE (in ppm) 3.50 3.23 



structure an anion. The force field (based on Am- 
ber) described previously was used for the MM 
region. DFT was used to describe the remaining 
102 atoms in domain 2. We used the LANL2DZ 
psuedopotential and basis sets as defined in Gaus- 
sian 09 with the BLYP functional. For computa- 
tional efficiency, we use density fitting with the 
W06 fitting basis 29J. It is evident that for this 
type of partition the gain in computational time, 
with respect to a full QM calculation, is minimal, 
since the QM region has effectively 24 gold atoms. 
Nevertheless, such partitioning scheme represents 
a proof of concept and will clearly become com- 
putationally more efficient for larger clusters and 
for cases where a reduced number of ligands are 
described quantum mechanically. 



2.2.2 Au 38 

As a second example, we use our QM/MM model 
to optimize the structure of Au38(SCH 3 ) 2 4. The 
X-ray structure was determined recently [30] along 
with a full DFT comparison 31. The Au 38 (SCH 3 )24 
structure has a bi-icosahedral Au2 3 core, consist- 
ing of two icosahedral Aui 3 units joined by three 
shared Au atoms at the center. This core is pro- 
tected by six long [(Au)2(SCH) 3 ] semi-ring (V- 
shape motifs) and three short [Au(SCH3)2] semi- 
rings (staple motifs) (see figure[3]). The entire struc- 
ture has D 3 symmetry, the largest Abelian point 



group being C2. 

By analogy with our partitioning of the [AU25 
(SCH 3 )i8]~ system, we describe the bi-icosahedral 
Au2 3 core in Au 3 g(SCH 3 ) 2 4 by our modified Am- 
ber force field. The outer shell comprising the V- 
shape and staple motifs was described quantum 
mechanically using DFT. As before, we use the 
BLYP pure functional with LANL2DZ and density 
fitting basis set W06. Charge of the QM region was 
assumed to be zero. 



2.2.3 Semi-rings 

We show in our final example, that one can use the 
same prescription proposed earlier to isolate a por- 
tion of the monolayer for a QM description, while 
retaining an MM description of the remainder. For 
this purpose, we took the [Au 2 5(SCH 3 )i 8 ]~ struc- 
ture and changed the ligands in two of the semi- 
rings from methyl to ethylphenyl, thus forming 
[Au25(SCH 3 ) 12 (SCH2CH2Ph) 6 ]- (see figure^ . The 
QM region contains two semi-rings on opposite 
sides of the structure and the symmetry of the 
original structure is preserved. Again the Aui 3 core 
and remaining (Au)2(SCH 3 ) 3 semi-rings are in the 
MM region. We optimized the geometry and calcu- 
lated NMR properties of 13 C and 1 H. Results are 
reported in the following section. 
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Table 5 Structure validation for different QM/MM models. All distances are reported in Angstroms. 



Cluster d,2i {d,2i)xRAY d rms (d rm s)x ray Ss (%) 

[Au2 5 (SCH 3 )i8]" 2.84 2.87 2.73 2.76 1.09 

[Au3 8 (SCH 3 )24] 2.09 2.02 1.71 1.76 2.84 

[Au25(SCH3)i2(SCH2CH 2 Ph)6]- 2.87 2.87 2.72 2.76 1.45 




Fig. 4 Hybrid structure [Au2 5 (SCH 3 )i2(SCH 2 CH2 
Ph)g] - . The two semi-rings, belonging to the QM re- 
gion, have ethyl-phenyls as ligands (shown by sticks). 
The remaining four semi-rings have methyls, which to- 
gether with the AU13 core belong to the MM region 
(shown in balls & sticks). 



3 Results 

We report in Table [5] an analysis of the optimized 
structure for the different clusters described in the 
previous section. To assess the accuracy of our 
structure, we define the parameters d 2 i and d rms 
as the distance between two neighboring Au atoms 
in the core and the "root mean square" distance 
of all Au atoms from the central atom in the core, 
respectively. We define the percent error 6$ as 

g _ \{d rms )xRAY — {drms) MODEL] ^ 
{drms) X BAY 

We note that the maximum error in structure 
optimization using our QM/MM models is ~ 3 %, 
which is remarkably good. We then used the op- 
timized structures to calculate NMR properties of 
13 C and 1 H, which we present in Tables [6] and [7J 
Guided by the benchmark studies of Cheeseman 
et al. 32 , our calculations of NMR chemical shifts 
are performed with the B3LYP hybrid functional, 



using the LANL2DZ basis set for Au, and the 6- 
31G(d,p) for S, C and H. We see that the NMR 
isotropic shieldings obtained using the QM/MM 
model compare very well to that found using full 
DFT on the same [Au 25 (SCH 3 )i 2 (SCH 2 CH 2 Ph) 6 ]- 
structure. We also calculate the mean unsigned er- 
ror (MUE) to compare theory with NMR experi- 
ments reported recently [10]. The mean unsigned 
error is defined as 

MUE — — ^(l^modei — ^experiment I) 
Til 

l 

where Uj is the total number of different types 
of C and H atoms, and 6 are the NMR chem- 
ical shifts. The maximum MUE is 3.5 ppm for 
13 C and 0.18 ppm for X H NMR. Thus we see that 
both 13 C and 1 H agree remarkably well with the 
known experimental results. Although small, we 
believe that the discrepancies arise from the fact 
that the experimental structure has six semi-rings 
with CH 2 CH 2 Ph as the functional group for all 
ligands whereas the QM/MM and full DFT mod- 
els have two semi-rings with CH 2 CH 2 Ph as ligands 
and the rest replaced with CH3 , so we are neglect- 
ing some of the electrostatic interactions that arise 
due to the difference in the structure. However, the 
error seems to be small in comparison to the huge 
gain in computational efficiency. More precisely, 
the time taken to caculate the NMR shielding ten- 
sors with the QM/MM model was approximately 1 
hour, whereas the same calculation with full DFT 
took 12 hours using the same hardware. 

4 Summary 

The use of X-ray crystallography to determine the 
structure of gold nanoparticles has been limited, 
owing partly to the difficulty of obtaining samples 
of sufficiently uniform size for the growth of sin- 
gle crystals. Until recently, there had been just one 
example of an Au cluster, Auio 2 (p-MBA)44, which 
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Table 7 1 H 


NMR chemical 


shifts 


S (in ppm) for 


[AU25(SCH 3 )12 


(SCH 2 CH 2 Pli) 6 


cluster. 


Atom 


QM/MM 


DFT 


H vriorirnDnf 1 M 
HjXjJfcJI llllclll _L U 


CH2 , inner 


2.92 


3.23 


q f\ h f ad 

0.0U 


OL— CH2, outer 


3.18 


2.97 


3.13 


CH2 ; inner 


2.93 


2.95 


3.13 


/3— CH2. outer 


3.03 


3.07 


t2.93 


111 CH inner 


7.28 


7.24 


7.15 


111 C^liouter 


7.29 


7.19 


7.19 


O^CHinner 


7.19 


7.54 


7.19 


_ CH ou ter 


7.23 


7.35 


7.14 


P~CH inner 


7.23 


7.02 


7.08 


P~CH G uter 


7.27 


7.09 


7.15 


MUE (in ppm) 0.18 


0.18 





had been determined by crystallography |6J. Since 
then, there have been successful crystal structure 
determinations of smaller sized gold clusters like 
AU25 [Z] and Au3g [30]. The recent availability of 
detailed structural data has created an opportu- 
nity for the construction and validation of com- 
putational methods for modeling these systems. 
Starting from the analysis of [Au 2 5(SCH 3 ) 18 ]~ struc- 
ture, we have developed generalized hybrid QM/MM 
models to accurately describe the structure of mono- 
layer protected gold clusters. We have found that 
errors in the calculated geometries from our QM/MM 
model are ~ 3%, compared to the crystal struc- 
tures. In addition, the isotropic chemical shifts ob- 
tained from NMR calculation agree very well among 
the QM/MM model, full DFT calculations and 
experiment. Use of these hybrid models not only 
saves computational resources but also enables the 
study of physical properties and chemical reactions 
for larger clusters where physical insight can be 
obtained by describing only a small region of the 
system at a high QM level. 
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